Lab 1: NumPy exercises¶

This ungraded assignment is composed of 26 NumPy exercises divided in three difficulty levels: beginner, intermediate, and advanced. In order to understand and solve the exercises, and more generally to gain the maximum benefit from this assignment, it is highly recommended that you read Chapter 1 in the course notes available on Blackboard.

IMPORTANT: Please refer to Chapter 1 as you work through this assignment.

Instructions¶

  • Download a copy of this notebook from Blackboard.

  • Run jupyter notebook and open the .ipynb file.

    • Keep the notebook inside the folder it was downloaded with. This folder contains all the dependencies needed for the notebook to work properly.
  • Work alone or with a partner to solve the quizzes.

    • You are supposed to fill in or modify the code marked with the comment # YOUR CODE HERE
    • You can check your answers by running the tests provided at the end of each quiz

Contents¶

  • Beginner level
Exercise Topic
Quiz 1.1 Array creation I
Quiz 1.2 Array creation II
Quiz 1.3 Attribute .shape
Quiz 1.4 Reshaping
Quiz 1.5 Vectorization I
Quiz 1.6 Vectorization II
Quiz 1.7 Keyword axis
Quiz 1.8 Simple indexing
Quiz 1.9 Slicing
  • Intermediate level
Exercise Topic
Quiz 2.1 Creation
Quiz 2.2 Slicing I
Quiz 2.3 Slicing II
Quiz 2.4 Reshaping I
Quiz 2.5 Reshaping II
Quiz 2.6 Reshaping III
Quiz 2.7 Indexing
Quiz 2.8 Broadcasting
Quiz 2.9 Scalar product
  • Advanced level
Exercise Topic
Quiz 3.1 Checkboard
Quiz 3.2 Column swapping
Quiz 3.3 Find minimum
Quiz 3.4 Sum pairs of cols
Quiz 3.5 Softmax
Quiz 3.6 Linear regression
Quiz 3.7 Pairwise distances
Quiz 3.8 Advanced indexing

Required packages¶

For this assignment, you need to import the following packages.

  • Numpy - The fundamental package for scientific computing with Python.
In [2]:
import numpy as np

Part 1. Beginner level¶

Quiz 1.1¶

Create a vector of length 5, filled with zeros.

Hint: np.zeros()

In [3]:
vector_zeros = np.zeros(5) # SOLUTION
In [ ]:
from numpy.testing import assert_almost_equal
assert_almost_equal(vector_zeros, [0, 0, 0, 0, 0])

Quiz 1.2¶

Create a matrix with 8 rows and 4 columns, filled with one.

Hint: np.ones()

In [7]:
matrix_ones = np.ones((8,4)) # SOLUTION
In [8]:
from numpy.testing import assert_almost_equal
assert_almost_equal(matrix_ones, 
[[1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1],
 [1, 1, 1, 1]])

Quiz 1.3¶

Get the number of rows and columns of a matrix.

Hint: .shape

In [9]:
# This is the array to work with
matrix = np.array([[ 1,  2,  3], 
                   [ 4,  5,  6], 
                   [ 7,  8,  9], 
                   [10, 11, 12]])

nrows = matrix.shape[0] # SOLUTION
ncols = matrix.shape[1] # SOLUTION
In [10]:
assert nrows == 4
assert ncols == 3

Quiz 1.4¶

Propose three different ways to reshape an array into a vector (1d array).

Hints:

  • .reshape()
  • .ravel()
  • .flatten()
In [16]:
# This is the array to work with
volume = np.array([[[5,6,3],
                    [3,2,7],
                    [3,5,1]],
                   
                   [[5,6,3],
                    [3,2,7],
                    [3,5,1]]])

flat_vector = volume.ravel() # SOLUTION
In [17]:
from numpy.testing import assert_almost_equal
assert_almost_equal(flat_vector, [5, 6, 3, 3, 2, 7, 3, 5, 1, 5, 6, 3, 3, 2, 7, 3, 5, 1])

Quiz 1.5¶

Given an array x, compute the following expression:

log(N−1∑i=0exp(xi)).log⁡(∑i=0N−1exp⁡(xi)).

Hint:

  • exp(⋅)exp⁡(⋅) denotes the exponential, which can be computed in numpy with the function np.exp().
  • ∑i⋅∑i⋅ denotes a summation, which can be computed in numpy with the function np.sum().
  • log(⋅)log⁡(⋅) denotes the logarithm, which can be computed in numpy with the function np.log().
In [19]:
# This is the array to work with
x = np.array([7.1, -5.7, 13])

log_sum_exp = np.log(np.exp(x).sum()) # SOLUTION
In [20]:
from numpy.testing import assert_almost_equal
assert_almost_equal(log_sum_exp, 13.00273570692086)

Quiz 1.6¶

Find the minimum and maximum values of a vector.

Hint:

  • .min()
  • .max()
In [21]:
# This is the vector to work with
vector = np.array([5, -1, 3, 12, 8, 0])

min_value = vector.min() # SOLUTION
max_value = vector.max() # SOLUTION

print("min:", min_value)
print("max:", max_value)
min: -1
max: 12
In [22]:
assert min_value == -1
assert max_value == 12

Quiz 1.7¶

Compute the mean value of a matrix, of its rows, and of its columns.

Hint:

  • .mean()
In [24]:
np.random.seed(1)

# This is the matrix to work with
matrix = np.random.randint(10,size=(3,4))

mean_all  = matrix.mean() # SOLUTION
mean_rows = matrix.mean(axis=1) # SOLUTION
mean_cols = matrix.mean(axis=0) # SOLUTION
In [32]:
from numpy.testing import assert_almost_equal
assert_almost_equal(mean_all, 4.66666666666666)
assert_almost_equal(mean_rows, [6.75, 2.0, 5.25])
assert_almost_equal(mean_cols, [3.66666667, 5.66666667, 4., 5.33333333])

Quiz 1.8¶

Swap (in-place) the 2nd and 4th element of a vector.

Remark: Single indexing extracts one element from an array. In Python, single values are always copied when assigned to a variable.

In [34]:
# This is the vector to work with
vector = np.array([5, -1, 3, 12, 8, 0])

print("--- before swap ---")
print(vector)

# YOUR CODE HERE
None
# BEGIN SOLUTION
vector[1], vector[3] = vector[3], vector[1]
# END SOLUTION
print("\n--- after swap ---")
print(vector)
--- before swap ---
[ 5 -1  3 12  8  0]

--- after swap ---
[ 5 12  3 -1  8  0]
In [36]:
from numpy.testing import assert_almost_equal
assert_almost_equal(vector, [5, 12, 3, -1, 8, 0])

Quiz 1.9¶

Given a matrix, extract the elements indicated in the figure below.

Hint: Recall that you can extract a contiguous part of a matrix via the slicing notation:

matrix[start_row : stop_row : step_row , start_col : stop_col : step_col]

You can omit any of the slicing indices. For example, these are valid slicing: start:stop, start:, :stop, :, ::step.

slicing.png

In [38]:
# This is the matrix to work with
matrix = np.array([[ 1,  2,  3], 
                   [ 4,  5,  6], 
                   [ 7,  8,  9], 
                   [10, 11, 12]])

A = matrix[0] # SOLUTION
B = matrix[:,0] # SOLUTION
C = matrix[2:,1:] # SOLUTION
D = matrix[::2] # SOLUTION
E = matrix[0:2,0:2] # SOLUTION
F = matrix[:,::2] # SOLUTION
G = matrix[1::2] # SOLUTION
H = matrix[:2,1:] # SOLUTION
In [46]:
from numpy.testing import assert_almost_equal
assert_almost_equal(A, [1, 2, 3])
assert_almost_equal(B, [ 1,  4,  7, 10])
assert_almost_equal(C, [[8, 9], [11, 12]])
assert_almost_equal(D, [[1, 2, 3], [ 7,  8,  9]])
assert_almost_equal(E, [[1, 2], [4, 5]])
assert_almost_equal(F, [[ 1,  3], [ 4,  6], [ 7,  9], [10, 12]])
assert_almost_equal(G, [[ 4,  5,  6], [10, 11, 12]])
assert_almost_equal(H, [[2, 3], [5, 6]])

Part 2. Intermediate level¶

Quiz 2.1¶

Create a vector of size 10 filled with zero, and set the first and last elements to 1.

Remark: The solution can take more than one line of code.

In [48]:
# ADD CODE HERE
head_tail = None
# BEGIN SOLUTION
head_tail = np.zeros(10)
head_tail[0] = 1
head_tail[-1] = 1
# END SOLUTION
print(head_tail)
[1. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
In [50]:
from numpy.testing import assert_almost_equal
assert_almost_equal(head_tail, [1, 0, 0, 0, 0, 0, 0, 0, 0, 1])

Quiz 2.2¶

Create a vector of size 20 with an alternating pattern of 0 and 1.

Expected output: [0, 1, 0, 1, 0, 1, ..., 0, 1]

In [51]:
# ADD CODE HERE
zero_one = None
# BEGIN SOLUTION
zero_one = np.zeros(20)
zero_one[1::2] = 1
# END SOLUTION
print(zero_one)
[0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
In [52]:
from numpy.testing import assert_almost_equal
assert_almost_equal(zero_one, [0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1, 0,1])

Quiz 2.3¶

Swap the first half of a vector with its last half.

Example: [1,2,3, 4,5,6] --> [4,5,6, 1,2,3]

Hints:

  • The size of a vector is given by the attribute .size or .shape[0].
  • Use the integer division n//2 to compute the length of half a vector.
  • The function np.concatenate() merges a copy of the input vectors into a new array.
In [54]:
# This is the vector to work with
vector = np.array([4,2,8,  7,1,3])

# ADD CODE HERE
swapped = None
# BEGIN SOLUTION
n = vector.size
head = vector[:n//2]
tail = vector[n//2:]
swapped = np.concatenate([tail,head])
# END SOLUTION
print("vector:", vector)
print("swapped:", swapped)
vector: [4 2 8 7 1 3]
swapped: [7 1 3 4 2 8]
In [56]:
from numpy.testing import assert_almost_equal
assert_almost_equal(vector,  [4, 2, 8, 7, 1, 3])
assert_almost_equal(swapped, [7, 1, 3, 4, 2, 8])

Quiz 2.4¶

Reshape an array of n elements into a matrix with n/2 rows and 2 columns.

Hint: .reshape()

In [57]:
# This is the array to work with
a = np.array([[5,6,3], [3,2,7], [4,2,8], [7,1,3]])

reshaped = a.reshape(a.size//2, 2) # SOLUTION

print("--- before ---")
print(a)
print("--- after ---")
print(reshaped)
--- before ---
[[5 6 3]
 [3 2 7]
 [4 2 8]
 [7 1 3]]
--- after ---
[[5 6]
 [3 3]
 [2 7]
 [4 2]
 [8 7]
 [1 3]]
In [59]:
from numpy.testing import assert_almost_equal
assert_almost_equal(reshaped, [[5, 6], [3, 3], [2, 7], [4, 2], [8, 7], [1, 3]])

Quiz 2.5¶

Given a vector with an even number of elements, compute the sum of two consecutive elements.

  • Expected result: Vector holding the sum of 1st and 2nd elements, then the sum of 3rd and 4th elements, and so on.
  • Example: [1, 2, 3, 4, 5, 6] --> [3, 7, 11]

Hint: Reshape the vector into a matrix with n/2 rows and 2 columns.

In [60]:
np.random.seed(6)

# This is the vector to work with
a = np.random.randint(10,size=(6,))

sum_pairs = a.reshape(a.size//2, 2).sum(axis=1) # SOLUTION

print("vector:", a)
print("sum by pairs:", sum_pairs)
vector: [9 3 4 0 9 1]
sum by pairs: [12  4 10]
In [62]:
from numpy.testing import assert_almost_equal
assert_almost_equal(sum_pairs, [12,  4, 10])

Quiz 2.6¶

Given a matrix with an even number of rows, compute the sum of two consecutive rows.

  • Expected result: Vector holding the sum of elements on 1st and 2nd rows, the sum of elements on 3rd and 4th rows, and so on.

Hint: Reshape the matrix so as to align the rows that must be summed together.

In [63]:
np.random.seed(1)

# This is the matrix to work with
A = np.random.randint(10,size=(6,4))

sum_rows = A.reshape(-1,2*A.shape[1]).sum(axis=1) # SOLUTION

print("---- A ----")
print(A)
print()
print("Sum:", sum_rows)
---- A ----
[[5 8 9 5]
 [0 0 1 7]
 [6 9 2 4]
 [5 2 4 2]
 [4 7 7 9]
 [1 7 0 6]]

Sum: [35 34 41]
In [65]:
from numpy.testing import assert_almost_equal
assert_almost_equal(sum_rows, [35, 34, 41])

Quiz 2.7¶

Replace the minimum value of a matrix by -10.

Hint: You can reshape the matrix into a vector, get the position of the minimum element with .argmin(), and modify the minimum value.

In [66]:
np.random.seed(1)

# This is the matrix to work with
matrix = np.random.randint(50, size=(4,6))

print("--- before ---")
print(matrix)

# ADD CODE HERE
None
# BEGIN SOLUTION
i = matrix.argmin()
matrix.flat[i] = -10
# END SOLUTION
print()
print("----- after -----")
print(matrix)
--- before ---
[[37 43 12  8  9 11]
 [ 5 15  0 16  1 12]
 [ 7 45  6 25 20 37]
 [18 20 11 42 28 29]]

----- after -----
[[ 37  43  12   8   9  11]
 [  5  15 -10  16   1  12]
 [  7  45   6  25  20  37]
 [ 18  20  11  42  28  29]]
In [68]:
from numpy.testing import assert_almost_equal
assert_almost_equal(matrix,
[[ 37,  43,  12,   8,   9,  11],
 [  5,  15, -10,  16,   1,  12],
 [  7,  45,   6,  25,  20,  37],
 [ 18,  20,  11,  42,  28,  29]])

Quiz 2.8¶

Divide each row of a matrix by its maximum value, and put the result in a new matrix. Exploit broadcasting to avoid loops!

Hint: Remember what the argument keepdims does in the function np.max().

In [69]:
np.random.seed(1)

# This is the matrix to work with
matrix = np.random.randint(10, size=(3,6))

# STEP 1. Compute the max along the rows
max_values = matrix.max(axis=1, keepdims=True) # SOLUTION

# STEP 2. Divide the matrix by the max values
new_matrix = matrix / max_values # SOLUTION
In [74]:
from numpy.testing import assert_almost_equal
assert_almost_equal(max_values, [[9], [9], [7]])
assert_almost_equal(new_matrix, 
[[0.55555556, 0.88888889, 1.        , 0.55555556, 0.        , 0.        ],
 [0.11111111, 0.77777778, 0.66666667, 1.        , 0.22222222, 0.44444444],
 [0.71428571, 0.28571429, 0.57142857, 0.28571429, 0.57142857, 1.        ]])

Quiz 2.9¶

Compute the following expression:

sigmoid(w⊤x)=11+e−(w0+w1x1+⋯+wQxQ)sigmoid(w⊤x)=11+e−(w0+w1x1+⋯+wQxQ)

where

w=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣w0w1⋮wQ⎤⎥ ⎥ ⎥ ⎥ ⎥⎦∈RQ+1x=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣1x1⋮xQ⎤⎥ ⎥ ⎥ ⎥ ⎥⎦∈RQw=[w0w1⋮wQ]∈RQ+1x=[1x1⋮xQ]∈RQ

with the convention x0=1x0=1.

Hint: Note that the vectors don't have the same size:

  • w - vector of shape (Q+1,)
  • x - vector of shape (Q,).

To get around this, you can proceed as follows:

  • Slice w so as to extract the elements w[1], ..., w[Q].
  • Compute the scalar product between x and the above elements.
  • Add w[0] to the final result.
In [75]:
# These are the vectors to work with
x = np.array([3, 2, 7])
w = np.array([-3, 2, 0.5, -1])

# STEP 1. Compute 'w0 + w1*x1 + ... + wQ*xQ'.
w_dot_x = w[0] + np.dot(x, w[1:]) # SOLUTION
    
# STEP 2. Compute the sigmoid of 'x_dot_w', as defined above.
logistic = 1/(1+np.exp(-w_dot_x)) # SOLUTION
In [77]:
from numpy.testing import assert_almost_equal
assert_almost_equal(logistic, 0.04742587317756678)

Part 3. Advanced level (optional)¶

Quiz 3.1¶

Create a 8x8 matrix and fill it with a checkerboard pattern.

Expected output:

[[0 1 0 1 0 1 0 1] [1 0 1 0 1 0 1 0] [0 1 0 1 0 1 0 1] [1 0 1 0 1 0 1 0] [0 1 0 1 0 1 0 1] [1 0 1 0 1 0 1 0] [0 1 0 1 0 1 0 1] [1 0 1 0 1 0 1 0]]

In [78]:
# ADD CODE HERE
checkboard = None
# BEGIN SOLUTION
checkboard = np.zeros((8,8))
checkboard[1::2,::2] = 1
checkboard[::2,1::2] = 1
# END SOLUTION
checkboard
Out[78]:
array([[0., 1., 0., 1., 0., 1., 0., 1.],
       [1., 0., 1., 0., 1., 0., 1., 0.],
       [0., 1., 0., 1., 0., 1., 0., 1.],
       [1., 0., 1., 0., 1., 0., 1., 0.],
       [0., 1., 0., 1., 0., 1., 0., 1.],
       [1., 0., 1., 0., 1., 0., 1., 0.],
       [0., 1., 0., 1., 0., 1., 0., 1.],
       [1., 0., 1., 0., 1., 0., 1., 0.]])
In [79]:
from numpy.testing import assert_almost_equal
assert_almost_equal(checkboard,
[[0, 1, 0, 1, 0, 1, 0, 1],
 [1, 0, 1, 0, 1, 0, 1, 0],
 [0, 1, 0, 1, 0, 1, 0, 1],
 [1, 0, 1, 0, 1, 0, 1, 0],
 [0, 1, 0, 1, 0, 1, 0, 1],
 [1, 0, 1, 0, 1, 0, 1, 0],
 [0, 1, 0, 1, 0, 1, 0, 1],
 [1, 0, 1, 0, 1, 0, 1, 0]])

Quiz 3.2¶

Swap (in-place) the 2nd and 4th column of a matrix.

Recall when data is shared or copied:

  • No copy - A simple assignment produces an alias of the original array.
  • Shared data - Slicing creates a view of the original array.
  • Copied data - The method .copy() makes a copy of the original array.
In [80]:
# This is the matrix to work with
matrix = np.arange(15).reshape(3,5)

print("--- matrix (before swap) ---")
print(matrix)

# ADD CODE HERE
# BEGIN SOLUTION
temp = matrix[:,1].copy()
matrix[:,1] = matrix[:,3]
matrix[:,3] = temp
# END SOLUTION

print("\n--- matrix (after swap) ---")
print(matrix)
--- matrix (before swap) ---
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]]

--- matrix (after swap) ---
[[ 0  3  2  1  4]
 [ 5  8  7  6  9]
 [10 13 12 11 14]]
In [82]:
from numpy.testing import assert_almost_equal
assert_almost_equal(matrix,
[[ 0,  3,  2,  1,  4],
 [ 5,  8,  7,  6,  9],
 [10, 13, 12, 11, 14]])

Quiz 3.3¶

Find the element in a vector nearest to the value 1.

Hint: Given the vector x, compute the absolute values of x-1, and then find the index i of the minimum element.

  • np.abs()
  • np.argmin()

Example: step-by-step

  • Consider the vector x = [5, -1, 3, 1.2, 8, 0]. What is the element in x nearest to the value 1?
  • The operation x-1 results in the vector [4, -2, 2, 0.2, 7, -1].
  • The absolute value of x-1 is thus [4, 2, 2, 0.2, 7, 1].
  • What is the minimum element of this vector? Take notice of its position. It is the same as the element you are looking for!
In [83]:
# This is the vector to work with
vector = np.array([5, -1, 3, 1.2, 8, 0])

# ADD CODE HERE
nearest_value = None
# BEGIN SOLUTION
abs_diff = np.abs(vector-1)
i = abs_diff.argmin()
nearest_value = vector[i]
# END SOLUTION
print('nearest value:', nearest_value)
nearest value: 1.2
In [84]:
assert nearest_value == 1.2

Quiz 3.4¶

Given a matrix with an even number of columns, compute the sum of two consecutive columns.

  • Expected result: Vector holding the sum of elements on 1st and 2nd columns, the sum of elements on 3rd and 4th columns, and so on.

Hint: Recall that a matrix is stored in row-major order, so you need to come up with a creative way of reshaping the matrix.

In [85]:
np.random.seed(1)

# This is the matrix to work with
A = np.random.randint(10,size=(6,4))

sum_cols = A.sum(axis=0).reshape(-1,2).sum(axis=1) # SOLUTION

print("--- A ---")
print(A)
print()
print("Sum:", sum_cols)
--- A ---
[[5 8 9 5]
 [0 0 1 7]
 [6 9 2 4]
 [5 2 4 2]
 [4 7 7 9]
 [1 7 0 6]]

Sum: [54 56]
In [88]:
from numpy.testing import assert_almost_equal
assert_almost_equal(sum_cols, [54, 56])

Quiz 3.5¶

  • axis=None - sum over all the elements
  • axis=0 - sum over the columns
  • axis=1 - sum over the rows.

Implement a function that computes the following expression:

softmax(x)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ex0∑N−1j=0exj⋮exN−1∑N−1j=0exj⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.softmax(x)=[ex0∑j=0N−1exj⋮exN−1∑j=0N−1exj].

When the input is a matrix, the function must operate according to an input parameter axis:

Hint: In the skeleton function given below, print the shapes of intermediate variables x_exp, x_sum and s to help you with broadcasting.

  • The shape of x_sum must be (1,1), (1,3) or (2,1) depending on the input parameter axis.

  • The shape of x_exp and s must always be (2,3).

  • Then, the division x_exp/x_sum will work due to broadcasting.

In [89]:
def softmax(x, axis=None):
    """Compute the softmax of x"""
        
    # Compute the exponential.
    x_exp = np.exp(x) # SOLUTION

    # Sum the elements along the specified axis. Don't forget to keep the dimensions!
    x_sum = np.sum(x_exp, axis=axis, keepdims=True) # SOLUTION
    
    # Compute the division. It should automatically use broadcasting!
    s = x_exp / x_sum # SOLUTION

    return s
In [92]:
x = np.array([[9, 2, 5],
              [7, 5, 0]])

s_all  = softmax(x)
s_cols = softmax(x, axis=0)
s_rows = softmax(x, axis=1)
In [98]:
from numpy.testing import assert_almost_equal

assert_almost_equal(s_all, 
[[8.52513572e-01, 7.77391752e-04, 1.56143307e-02],
 [1.15375166e-01, 1.56143307e-02, 1.05208533e-04]])

assert_almost_equal(s_cols,
[[0.88079708, 0.04742587, 0.99330715],
 [0.11920292, 0.95257413, 0.00669285]])

assert_almost_equal(s_rows,
[[9.81135202e-01, 8.94679497e-04, 1.79701181e-02],
 [8.80090205e-01, 1.19107257e-01, 8.02538386e-04]])

Quiz 3.6¶

Implement a function that computes the following expression:

J(w)=N∑n=1(y(n)−w⊤x(n))2,J(w)=∑n=1N(y(n)−w⊤x(n))2,

where

w=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣w0w1⋮wQ⎤⎥ ⎥ ⎥ ⎥ ⎥⎦∈RQ+1x(n)=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣1x1⋮xQ⎤⎥ ⎥ ⎥ ⎥ ⎥⎦∈RQy(n)∈Rw=[w0w1⋮wQ]∈RQ+1x(n)=[1x1⋮xQ]∈RQy(n)∈R

with the convention x(n)0=1x0(n)=1.

Hint: This computation can be vectorized as follows.

  • The vectors x(n)x(n) are stacked in a matrix X∈RN×QX∈RN×Q, whereas the scalars y(n)y(n) are stacked in a vector y∈RNy∈RN:

X=⎡⎢ ⎢ ⎢⎣__x(1)⊤__⋮__x(N)⊤__⎤⎥ ⎥ ⎥⎦y=⎡⎢ ⎢ ⎢⎣y(1)⊤⋮y(N)⊤⎤⎥ ⎥ ⎥⎦.X=[__x(1)⊤__⋮__x(N)⊤__]y=[y(1)⊤⋮y(N)⊤].

  • The products z(n)=w⊤x(n)z(n)=w⊤x(n), for every index nn, are computed by multiplying XX by ww:

z=Xw=⎡⎢ ⎢⎣w⊤x(1)⋮w⊤x(N)⎤⎥ ⎥⎦.z=Xw=[w⊤x(1)⋮w⊤x(N)].

  • The distance JJ is now equal to the squared norm of the difference between zz and yy:

J=∥z−y∥2=N∑n=1(z(n)−y(n))2.J=‖z−y‖2=∑n=1N(z(n)−y(n))2.

Note hovewer that the shapes of X and w don't align for a matrix-vector product. To get around this, you can proceed as follows:

  • Slice w so as to extract the elements w[1], ..., w[Q].
  • Compute the product between X and the above elements.
  • Add w[0] to the final result.
In [99]:
def ssd_cost(w, X, y):
    """Calculates the sum of squared differences.
     
    Arguments:
    w -- flat vector of shape (Q+1,)
    X -- matrix of shape (N, Q)
    y -- flat vector of shape (N,)

    Returns:
    s -- Scalar
    """
    
    # Compute w0 + w1*X1 + ... + wQ*XQ
    w_dot_x = w[0] + np.dot(X, w[1:]) # SOLUTION

    # Compute the sum of squared distances between x_dot_w and y
    J = np.sum( (y - w_dot_x)**2 ) # SOLUTION
    
    return J
In [101]:
w = np.array([-3, 2, 0.5, -1])
X = np.array([[4, -5,  2],
              [5,  3, -3],
              [5,  2,  1],
              [8,  6,  9],
              [5, -6,  2]])
y = np.array([6, -3, 0.5, 1, 0.])

assert ssd_cost(w, X, y) == 322.75

Quiz 3.7¶

Compute the pairwise distance between the rows of two matrices.

X=⎡⎢ ⎢ ⎢ ⎢⎣__x1____x2__⋮__xN__⎤⎥ ⎥ ⎥ ⎥⎦Y=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣__y1____y2__⋮__yK__⎤⎥ ⎥ ⎥ ⎥ ⎥⎦dist(X,Y)=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣||x1−y1||||x1−y2||…||x1−yK||||x2−y1||||x2−y2||…||x2−yK||⋮⋮⋱⋮||xN−y1||||xN−y2||…||xN−yK||⎤⎥ ⎥ ⎥ ⎥ ⎥⎦X=[__x1____x2__⋮__xN__]Y=[__y1____y2__⋮__yK__]dist⁡(X,Y)=[||x1−y1||||x1−y2||…||x1−yK||||x2−y1||||x2−y2||…||x2−yK||⋮⋮⋱⋮||xN−y1||||xN−y2||…||xN−yK||]

Hint: It is possible to compute the differences X[i] - Y[j] via broadcasting. Here is how to do so.

  • Reshape the matrix X to size (N,1,M).

  • Reshape the matrix Y to size (1,K,M).

  • The reshaped arrays now align for broadcasting.

Remember: You can use the None operator to insert an axis of length 1.

  • The resulting array D is such that D[i,j] holds the vector X[i] - Y[j].

  • You can use the function np.linalg.norm(..., axis=...) to compute the norm along a given axis.

pairwise.png

In [102]:
def distance_matrix(X, Y):
    """Compute the distance between X[i] and Y[j] for all i,j"""
    
    # Reshape 'X' to the appropriate shape
    X_3D = X[:,None] # SOLUTION
    
    # Reshape 'Y' to the appropriate shape
    Y_3D = Y[None,:] # SOLUTION

    # Subtract the matrices. It should automatically use broadcasting.
    diff = X_3D - Y_3D # SOLUTION
    
    # Compute the norm  along the right axis. You should get a matrix of shape (X.shape[0], Y.shape[0])
    dist = np.linalg.norm(diff, axis=2) # SOLUTION

    return dist
In [104]:
np.random.seed(24)

X = np.random.randint(20,size=(4,5))
Y = np.random.randint(20,size=(3,5))

dist = distance_matrix(X, Y)

from numpy.testing import assert_almost_equal
assert_almost_equal(dist,
[[19.54482029, 11.48912529, 24.85960579],
 [23.87467277, 21.21320344, 18.92088793],
 [23.47338919, 24.75883681, 27.03701167],
 [18.24828759, 16.21727474, 12.36931688]])

Quiz 3.8¶

The Game of Life is a grid of square cells that can be in one of two possible states: live or dead.

life.png

The player creates an initial configuration of live cells and observes how they evolve. Every cell interacts with its 8 neighbours (the cells horizontally, vertically, or diagonally adjacent). At each step in time, the following transitions occur.

  • Underpopulation: Any live cell with fewer than two live neighbours dies.

  • Overcrowding: Any live cell with more than three live neighbours dies.

  • Stability: Any live cell with two or three live neighbours lives unchanged to the next generation.

  • Rebirth: Any dead cell with exactly three live neighbours becomes a live cell.

The initial pattern constitutes the 'seed' of the system. The first generation is created by applying the above rules simultaneously to every cell in the seed. Births and deaths happen simultaneously. The rules continue to be applied repeatedly to create further generations.

The Game of Life can be implemented by a matrix whose elements represent the single cells:

  • 0 means the cell is dead,
  • 1 means the cell is alive.
board = [[0,0,0,0,0,0],
         [0,0,0,1,0,0],
         [0,1,0,1,0,0],
         [0,0,1,1,0,0],
         [0,0,0,0,0,0],
         [0,0,0,0,0,0]]

The cells evolve through a for-loop which involves two operations at each iteration:

  • counting the neighbors,
  • applying the evolution rules.

Here is how the main function looks like.

def game_of_life(board, steps):
    for i in range(steps):
        count = neighbor_counting(board)
        board = evolution_rules(board, count)    
    return board

We provide you with the function neighbor_counting. Your job is to implement evolution_rules.

In [105]:
def neighbor_counting(board):
    padded = np.pad(board, (1,1), 'constant')
    count  = padded[ :-2, :-2] + padded[:-2, 1:-1] + padded[ :-2, 2:] \
           + padded[1:-1, :-2] +                   + padded[1:-1, 2:] \
           + padded[2:  , :-2] + padded[2: , 1:-1] + padded[2:  , 2:]
    return count

Implement a function that applies the evolution rules of the Game of Life (see the description above).

Hint: Use advanced indexing to avoid loops!

  • Translate the evolution rules to boolean conditions,
  • Use the boolean masks to set the new values.
In [106]:
def evolution_rules(board, count):
    """
    Arguments:
     board -- matrix whose elements indicate whether a cell is dead (0) or alive (1)
     count -- matrix whose elements indicate the number of alive neighbors of a cell

    Returns:
     new_board -- matrix with the same shape as 'board'
    """
    
    new_board = np.zeros_like(board)
            
    # Underpopulation
    R1 = (board == 1) & (count < 2) # SOLUTION 
    new_board[R1] = 0 # SOLUTION
    
    # overcrowding
    R2 = (board == 1) & (count > 3) # SOLUTION 
    new_board[R2] = 0 # SOLUTION
    
    # stability 
    R3 = (board == 1) &((count == 2) | (count == 3)) # SOLUTION 
    new_board[R3] = board[R3] # SOLUTION
    
    # rebirth
    R4 = (board == 0) & (count == 3) # SOLUTION 
    new_board[R4] = 1 # SOLUTION

    return new_board
In [107]:
board = np.array([[0,0,0,0,0,0],
                  [0,0,0,1,0,0],
                  [0,1,0,1,0,0],
                  [0,0,1,1,0,0],
                  [0,0,0,0,0,0],
                  [0,0,0,0,0,0]])

count = neighbor_counting(board)

new_board = evolution_rules(board, count)

from numpy.testing import assert_almost_equal
assert_almost_equal(new_board,
[[0, 0, 0, 0, 0, 0],
 [0, 0, 1, 0, 0, 0],
 [0, 0, 0, 1, 1, 0],
 [0, 0, 1, 1, 0, 0],
 [0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0]])

Bonus :-)¶

Congratulations! You made it to the end. As a reward, sit back and watch The Game of Life.

Run the following cells to play the animation. The creation process may take a while.

In [108]:
def game_of_life(board, steps):
    history = [board]
    for i in range(steps):
        count = neighbor_counting(board)
        board = evolution_rules(board, count)   
        history.append(board)
    return history
In [115]:
def create_animation(history, speed=1):
    import matplotlib.pyplot as plt
    import matplotlib.animation as animation

    fig = plt.figure(figsize=(7,7))
    im = plt.imshow(history[0], cmap=plt.cm.seismic, interpolation='bicubic')
    plt.axis('off')
    plt.close(fig)

    def animate(i):
        im.set_data(history[i])
        #im.autoscale()
        return (im,)

    anim = animation.FuncAnimation(fig, animate, frames=range(0,len(history),speed), interval=100, blit=True)
    return anim
In [116]:
glider_gun =\
[[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
 [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0],
 [0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1],
 [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1],
 [1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [1,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0],
 [0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
 [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
 [0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]]

seed = np.zeros((50, 70))
seed[1:10,1:37] = glider_gun
In [117]:
history = game_of_life(seed, steps=500)
In [118]:
from IPython.display import HTML

animation = create_animation(history, speed=2)

HTML(animation.to_jshtml())
Out[118]:
No description has been provided for this image

Reaction and diffusion of chemical species can produce a variety of patterns, reminiscent of those often seen in nature. The code below implements a more advanced version of The Game of Life, based on the Gray-Scott equations. Just execute the cell, no further action is required from you.

In [76]:
def create_seed(n):
    U = np.zeros((n+2, n+2))
    V = np.zeros((n+2, n+2))
    r = n//5 #20
    U[1:-1,1:-1] = 1.0
    U[n//2-r:n//2+r, n//2-r:n//2+r] = 0.50
    V[n//2-r:n//2+r, n//2-r:n//2+r] = 0.25
    return U,V


def laplacian(U, V):
    Lu = (                  U[0:-2, 1:-1] +
          U[1:-1, 0:-2] - 4*U[1:-1, 1:-1] + U[1:-1, 2:] +
                            U[2:  , 1:-1])
    Lv = (                  V[0:-2, 1:-1] +
          V[1:-1, 0:-2] - 4*V[1:-1, 1:-1] + V[1:-1, 2:] +
                            V[2:  , 1:-1])
    return Lu, Lv


def reaction(U, V, Lu, Lv, Du, Dv, F, k):
    u = U[1:-1, 1:-1]
    v = V[1:-1, 1:-1]
    uvv = u*v*v
    UU = np.zeros(U.shape)
    VV = np.zeros(V.shape)
    UU[1:-1, 1:-1] = u + (Du*Lu - uvv +  F   *(1-u))
    VV[1:-1, 1:-1] = v + (Dv*Lv + uvv - (F+k)*v    )
    return UU, VV
    
    
def chemical_reaction(n, steps, mode='bacteria-1'):
    
    if mode == 'bacteria-1':
        Du, Dv, F, k = 0.16, 0.08, 0.035, 0.065
    elif mode == 'bacteria-2':
        Du, Dv, F, k = 0.14, 0.06, 0.035, 0.065
    elif mode == 'coral':
        Du, Dv, F, k = 0.16, 0.08, 0.060, 0.062
    elif mode == 'fingerprint':
        Du, Dv, F, k = 0.19, 0.05, 0.060, 0.062
    elif mode == 'spirals':
        Du, Dv, F, k = 0.10, 0.10, 0.018, 0.050
    elif mode == 'spirals-dense':
        Du, Dv, F, k = 0.12, 0.08, 0.020, 0.050
    elif mode == 'spirals-fast':
        Du, Dv, F, k = 0.10, 0.16, 0.020, 0.050
    elif mode == 'unstable':
        Du, Dv, F, k = 0.16, 0.08, 0.020, 0.055
    elif mode == 'worms-1':
        Du, Dv, F, k = 0.16, 0.08, 0.050, 0.065
    elif mode == 'worms-2':
        Du, Dv, F, k = 0.16, 0.08, 0.054, 0.063
    elif mode == 'zebrafish':
        Du, Dv, F, k = 0.16, 0.08, 0.035, 0.060
    else:
        raise ValueError('Mode not recognized!')
    
    # seeds
    U, V = create_seed(n)
    
    # bookkeeping
    history = [V]
    
    # iterate
    for i in range(steps):
        Lu, Lv = laplacian(U, V)
        U, V = reaction(U, V, Lu, Lv, Du, Dv, F, k)
        history.append(V)
        
    return history

The function chemical_reaction() allows you to simulate a chemical reaction, much like in The Game of Life.

In [77]:
history = chemical_reaction(200, steps=5000, mode='bacteria-2')

Execute the next cell and enjoy the amination!

  • There are different reaction types that you can choose from. Look inside the function chemical_reaction() for the complete list.

  • You can go back, change the parameter steps or mode, and run again the simulation.

In [78]:
animation = create_animation(history, speed=50)

HTML(animation.to_jshtml())
Out[78]:
No description has been provided for this image